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ABSTRACT:  In  this  paper  we  illustrate  how  to  solve 
the  general  Helmholtz  equation  starting  from 
Laplace's  equation.  The  interesting  point  is  that  the 
Helmholtz  equation  has  a frequency  term  whereas 
Laplace’s  equation  is  the  static  solution  of  the  same 
boundary  value  problem.  In  this  new  formulation  the 
frequency  dependence  is  manifested  in  the  form  of  an 
excitation.  A new  boundary  integral  method  for 
solving  the  general  Helmholtz  equation  is  developed. 
This  new  formulation  is  developed  for  the  two- 
dimensional  Helmholtz  equation  with  the  method  of 
moments  Laplacian  solution.  The  main  feature  of  this 
new  formulation  is  that  the  boundary  conditions  are 
satisfied  independent  of  the  region  node 
discretizations.  The  numerical  solution  of  the  present 
method  is  compared  with  finite  difference  and  finite 
element  solutions  of  the  same  problem.  Application  of 
this  method  is  also  presented  for  the  computation  of 
cut-off  frequencies  for  some  canonical  waveguide 
structures. 

L INTRODUCTION 

The  two  dimensional  Helmholtz’s  equation  is  an 
important  equation  to  be  solved  in  many  numerical 
electromagnetics  problems  such  as  waveguide  related 
problems  and  appears  in  a variety  of  physical 
phenomena  and  engineering  applications,  such  as, 
acoustic  radiation  [1],  heat  conduction  [2],  and  water 
wave  propagation  [3],  In  semiconductor  device 
modeling,  Helmholtz’s  equation  arises  frequently  as 
an  intermediate  step  in  the  solution  of  the  nonlinear 
Poisson’s  problem.  To  solve  these  problems  diverse 
numerical  methods  have  been  reported  which  include 
finite  difference  [4],  finite  element  [5],  and  boundary 
integral  methods  (BIM)  [6-8].  Using  these 


conventional  methods,  it  has  been  found  that  fine  grids 
and  a large  number  of  elements  must  be  employed  to 
get  satisfactory  accuracy  [3].  This  requires  large 
computer  core  storage,  and  more  computational  time 
especially  for  the  iteration  scheme  of  the  nonlinear 
Poisson’s  problem  where  the  value  at  each  grid  point 
needs  to  be  updated  at  each  step  of  the  iteration. 
Further,  the  BIM  formulations  are  in  most  cases 
limited  to  homogeneous  Helmholtz  equation  and  tied 
closely  to  the  particular  problem  at  hand  [6].  In  this 
paper,  a simple  approach  to  solve  the  homogeneous 
and  nonhomogeneous  Helmholtz  equations  is 

proposed.  The  technique  is  based  on  the  computation 
of  Laplacian  potential  by  the  method  of  moments 
(MoM)  [9],  without  resorting  to  different  formulations 
using  Hankel  functions  as  it  is  commonly  done  in  BIM 
[10].  Besides  its  generality  to  solve  Laplace’s, 
Poisson’s,  and  Helmholtz’s  equations  in  one  single 
code  implementation,  the  present  method  will 
considerably  reduce  the  number  of  domain  grids 
compared  to  the  finite  difference  methods  and  does 
not  require  any  interpolation.  The  accuracy  of  the 
MOM  solution  from  this  new  formulation  will  be 
compared  with  the  solutions  of  fine  difference  method 
(FDM)  and  finite  element  method  (FEM),  using  the 
ELLPACK  implementation  [4]. 

II.  MATHEMATICAL  FORMULATION 

Consider  the  following  two-dimensional  elliptic 
equation  for  a smooth  function  defined  in  a 2-D 
region  defined  by  91  which  is  bounded  by  a contour  C 
so  that 

V24/U,>;)  + >.(x,y)'P(^,};)=FU,y)  (1) 
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where  X and  F are  known  functions  in  the  domain  9L 
The  general  form  of  (1)  includes,  as  specializations, 
the  following  cases: 

1)  Laplace’s  equation,  with  X = 0 and  F = 0 

2)  Poisson’s  equation,  with  X = 0 and  F * 0. 

3)  Helmholtz’s  equation,  with  X * 0 and  F * 0. 
So  only  one  general  formulation  addresses  the  solution 
to  all  of  the  three  cases  above.  In  addition  we  illustrate 
how  to  use  the  frequency  independent  solution  of  the 
Laplace’s  equation  to  solve  the  general  Helmholtz 
equation.  On  the  contour  C the  boundary  condition  can 
be  of  Dirichlet,  Neumann,  or  mixed  type,  as  given  by 
the  general  form 


(x'F+p— =y  (2) 

on 

where  a,  (3,  y,  are  known  spatial  functions  and 
represents  the  normal  derivative.  It  is  implied  that  a 
consistent  set  of  boundary  conditions  has  been  chosen 
for  the  problem. 

The  proposed  scheme  to  solve  the  given  boundary 
value  problem  starts  by  assuming  the  term  X 'F  to  be  a 
known  function,  and  including  it  with  the  given 
excitation  function  F,  and  thereby  reducing  (1)  to  the 
familiar  Poisson’s  equation 

v2'T(^y)=-©U,y)  (3) 

with 

0(x,;y)  = /l(*,;y)¥(x,;y)-F(*,;y).  (4) 


<t>Sx>  y) 


— {[©(*',  yVn 

2.71  „ 


K 


yj(x-x)2  +(j-/)2 


px'dy' 


(8) 


where  the  spatial  sets  ( x , y ) and  (x\  y')  denote  the 
spatial  coordinates  of  the  field  and  source  coordinates, 
respectively,  and  K is  an  arbitrary  constant.  Here, 
£n(  K ) represents  the  value  of  the  scalar  potential 
at  infinity  for  the  two  dimensional  case.  For  the  3D 
case,  this  term  does  not  exist,  as  the  potential  at 
infinity  is  zero.  Hence  the  Green’s  function  is  simply 
1/R  instead  of  the  natural  log  function.  The  value  of 
the  parameter  K is  chosen  to  be  100  for  the  2-D 
problem  as  the  reference  potential  at  infinity  is  not 
zero  but  finite. 

An  expression  similar  to  (8)  can  be  derived  to 
approximate  the  Laplacian  potential  <(>/,.  The  potential 
<})*  can  be  assumed  to  be  produced  by  some  equivalent 
sources  consisting  of  electrical  charges,  o located  on 
the  contour  C [12].  Then  the  potential  <()/,  at  any  point 
(, x , y)  can  be  obtained  from  o(x , y)  using  the  following 
integral  [9] 


i 


f 


— | a{x\y')l  n 

2k  c 


K 


V (-C— jc')2  +( y— J 


\dV 


(9) 


The  solution  to  Poisson’s  equation  in  (3)  can  be 
expressed  as 


(5) 

where  is  the  solution  to  the  homogeneous  Poisson’s 
equation  (Laplace’s  equation) 


where  /'  is  the  arc  length  on  the  contour  C. 

The  boundary  condition  of  the  homogeneous 
potential  <|)/f  is  obtained  from  (2)  and  (5)  as 


aA  + P 


dn 


a<t>e+P 


dn 


GO) 


o 

11 

CnJ 

> 

(6) 

and  is  the  particular  integral,  i.e., 

V = — 0(x,  y) 

(7) 

Here  we  use  the  particular  solution  of  the  Poison’s 
equation  given  by  [12] 


It  can  be  seen  that  (6)  along  with  the  boundary 
condition  of  (10),  constitute  a similar  boundary- value 
problem.  A similar  problem  was  addressed  in  [9] 
almost  thirty  years  ago. 

This  completes  the  formulation  of  the  problem. 
The  characteristic  features  of  this  formulation  are: 

a.  The  frequency  term  appears  in  an  explicit 
form. 

b.  The  Green’s  function  and  the  unknowns  are 
independent  of  frequency.  The  Green’s 
function  is  in  fact  the  static  Green’s 
function. 
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Next  we  illustrate  how  to  solve  these  coupled 
equations  both  for  the  homogeneous  part  and  the 
particular  integral. 

III.  NUMERICAL  SOLUTION  USING 
THE  METHOD  OF  MOMENTS  (MOM) 


c,(x,y) 


fen 

AC, 


K 


yl(x-x')2  + (y-y')2 


W 


(15) 


In  order  to  solve  the  above  equations  given 
by  (8)-(10)  we  employ  the  Method  of  Moments. 

To  evaluate  the  above  integral  in  (8)  we 
divide  the  domain  9t  into  N sub-regions.  The 
midpoint  coordinates  of  each  of  the  sub-regions  A9t, 
are  denoted  by  (. x y2i)-  Then  the  potential  §p  at  the 
field  point  (x,  y)  is  approximated  by 


Using  (4),  (11)  and  (14),  the  Helmholtz  potential  at 
an  arbitrary  field  point  ( x , y ) can  now  be  expressed  as 

^(x,y)  = 

M N 

Yi(Jici(x,y)+'£(At'¥i-Fi)ai(x,  y) 

i = 1 i = 1 

(16) 


<t>P  (x,  y)  = ]£  Q(x2i , y2i  )a,  (x,y)  (11) 


i = 1 


where 

ai(x,y)  ■ 


A9t, 


K 


j(x-x')2+(y-y')2 


dxdy 


(12) 


Hence  it  is  assumed  that  0(x,  y ) is  constant  within 
each  sub-region  A9t*  and  is  equal  to  the  value  0 

fe/j2/). 

Next  we  show  how  to  numerically  evaluate  (9). 
Here,  we  take  recourse  to  MOM.  We  expand  the 
surface  unknown  o using  a pulse  expansion.  For 
purposes  of  illustration  and  simplicity  let  us  choose 
point  matching.  Through  the  use  of  a pulse-expansion 
and  point-matching  techniques  we  will  solve  the 
present  problem  given  by  (9).  Furthermore,  if  the 
contour  C is  segmented  by  M straight  lines  of  length 
AC,  between  points  i and  i + 1,  then  a can  be 
represented  by  the  step  approximation 

= (13) 

i = 1 


where  for  simplicity  we  have  used  the  following 
abbreviations  in  (4); 

h = ^ (*2 h yn) 

F ,■  = F (x2/,  y2i) 

%=  ^ (*2i  y2d 

The  unknown  function  ¥(x,  y)  in  (16)  of  Helmholtz 
equation  can  now  be  evaluated  once  the  unknown 
terms,  c2  and  are  determined.  Next  a system  of 
two  matrix  equations  is  derived  and  solved  for  the 
unknowns  g;  and  *P/.  The  first  matrix  equation  of  this 
system  is  readily  obtained  by  satisfying  (16)  at  the 
midpoints  (x2>  yij)  of  each  of  the  N sub-region  A9^. 
Using  matrix  notation,  we  obtain 

l'*'- ] = [pp ]k]+ Ur, ] - Fi ] a?) 

where  p;i-  = cfay,  y2j),  and  qp  = afaq,  y2j). 

The  second  matrix  equation  is  now  obtained  by 
enforcing  the  boundary  condition  of  (10)  as  follows. 

We  define  fj  - ( x yj\  j = 1,  2,  ...,  M,  to  be  the 

midpoints  of  AC;-.  The  boundary  conditions  are 
enforced  at  each  fj . Substitution  of  (1 1)  and  (14)  into 

(10)  gives  the  following  set  of  equations 


where  P,(/)  is  a pulse  function  equal  to  1 on  AC/  and 
zero  elsewhere  and  a,  is  its  unknown  amplitude. 
Substituting  (13)  into  (9),  we  obtain  an  approximation 
for  <>A 

M 

</>h(x,y)  = YdOici(x,y)  (14) 

i = 1 

where 


M N 

4 -EM*, -*;)», 


i = 1 


i = 1 


where  7 j =7  ( rj  ) and 

an 


wji  = 


(18) 


(x,y)=fj 


(19) 
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( Q da  y 

bx  = a^+p-*-  (20) 

V on  )(XtyH 

Equation  (18)  can  now  be  written  in  matrix  notations 
as 

<21> 

Observe  that  (17)  and  (21)  form  a system  of  two 
equations  in  two  unknowns  o,  and  4 from  which 
they  can  be  solved.  We  use  (21)  to  obtain  an 
expression  for  a and  then  substitute  it  into  (17),  and 
after  simple  matrix  manipulations  we  obtain  the 
following  equation  for 

[a][>f,]=[b]  (22) 

where 

(23) 

-W 

and  [/]  denotes  the  V x V identity  matrix. 

From  the  last  three  equations  some  general 
comments  can  be  made: 

1)  Analytic  expressions  can  easily  be  derived 
for  the  evaluation  of  all  the  terms  of  the 
matrices  [A]  and  [B]  that  include  the  integrals 
(12)  and  (15).  This  is  true  at  least  for  the  2-D 
case. 

2)  The  matrix  elements  [w;7],  which  are  a part  of 
the  system  matrix  of  the  moment  method,  are 
similar  to  the  matrix  elements  obtained  in  [9]. 
For  different  Helmholtz  problems  with 
different  boundary  conditions,  [wjj]  remains 
unchanged.  This  observation  is  important  in 
an  iteration  scheme  where  the  boundary 
conditions  are  kept  the  same. 

3)  Problems  with  multiple  right-hand  sides  are 
solvable  with  minimum  additional 
computation  time  since  the  function  F in  (17) 
appears  only  as  a term  of  matrix  [5]  in  (24). 

4)  The  domain  and  the  contour  discretization 
schemes  are  totally  independent.  This  feature 
can  offer  the  flexibility  to  handle  boundary 
discontinuities  without  the  need  for  excessive 
domain  grid  generations.  This  will  be 


illustrated  through  numerical  examples  in  the 
next  section. 

5)  Once  ¥,■  is  determined  at  the  grid  points, 
using  Gaussian  elimination  for  instance  the 
potential  at  any  other  point  is  obtained  using 
ordinary  matrix  multiplications  as  shown  by 
(17).  No  interpolations  are  needed. 

Next,  we  illustrate  these  points  through  some 
numerical  examples. 

III.  STATIC  APPLICATIONS 

Example  l— Water  Wave  Propagation:  As  a first 
example,  we  consider  the  problem  of  water  wave 
propagation  in  a rectangular  basin  100m  x 100m  [3]. 
Denoting  by  ¥ the  water  evaluation,  then  the  wave 
propagation  is  governed  by  (1)  with  F = 0,  and  X = fc2 
where  k is  the  wave  number.  The  boundary  conditions 
used  are  displayed  in  Figure  1.  The  solution  by  MOM 
is  compared  to  the  solutions  obtained  by  finite 
difference  method  (FDM)  and  finite  element  method 
(FEM)  for  various  mesh  sizes  inside  the  domain.  The 
minimum  number  of  nodes  required  for  each  method 
to  converge  to  the  exact  solution  at  any  arbitrary  point 
is  optimized.  Figure  2 compares  the  solutions 
computed  along  the  line  x = 90  by  the  Method  of 
Moments  utilizing  (lxl),  (2x2)  and  (4x4)  grids. 
Figure  3 shows  the  solutions  obtained  by  the  finite 
element  method  utilizing  (5x5),  (10x10),  and  (20x20) 
grids.  Figure  4 shows  the  results  obtained  by  the  finite 
difference  method  utilizing  (5x5),  (15x15),  and 
(35x35)  grids.  In  all  these  figures,  the  solid  line 
represents  the  exact  solution.  Referring  to  these 
figures,  one  can  easily  observe  that  a considerably 
smaller  number  of  nodes  is  employed  with  the  present 
method  than  with  the  conventional  methods.  The 
minimum  number  of  the  domain  nodes  required  to 
obtain  a satisfactory  accuracy  using  MOM  for  the 
present  method  is  16  (4x4),  compared  to  400  (20x20) 
for  FEM  or  up  to  1000  nodes  for  FDM. 

Example  2 — MOST  Modeling : The  nonlinear 

Poisson’s  equation  plays  a key  role  in  numerical 
modeling  of  semi-conductor  devices.  Many  important 
characteristics  of  VLSI  devices  can  be  extracted  from 
the  solution  of  Poisson’s  equation.  The  most  common 
approach  to  the  numerical  solution  of  the  nonlinear 
Poisson’s  equation  is  based  on  the  application  of 
Newton’s  method  to  simultaneous  discretized 
equations  [15].  This  approach  often  requires  large 
storage  especially  for  fine  meshes,  as  is  the  case  for 
two-dimensional  modeling  of  the  MOST  [14]. 
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<t>  =f,(x) 


*T 
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y i 

k 

X 

► 

II 
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$=fb(x) 

cos  (ax)  if  20  < x < 80 
where  f (x)  = —2  if  X < 80 

1 if  x > 80 

fs(x)  = sin(ay),  fb(x)  = cos(ax),  and  a~2n 


Figure  1.  Physical  structure  and  the  associated  boundary  conditions 
associated  with  Example  1. 


In  this  section  the  MOM  is  applied  to  solve  the 
nonlinear  Poisson’s  equation  that  arises  in  the  MOST 
modeling.  We  consider  the  MOST  structure  of  Figure 
5 made  on  a p-type  substrate  with  doping  NA.  Under 
the  low  current  approximation  the  potential  O is 
governed  by  the  Poisson’s  equation  [16] 


Figure  2.  Solution  obtained  by  the  present  method. 


d2x¥  B2vF 

dx2  By2 


ft; 


n ; 


7 J 


(25) 


where  ¥ = ®/VT  is  the  normalized  potential,  VT  is  the 
thermal  potential,  n,  is  the  intrinsic  carrier 
concentration,  and  LD  is  the  Debye  length. 


$=V0-Vft 


o 

II 


$ =0 

Figure  5.  The  physical  structure  and  the  associated  boundary 
condition  for  MOST  modeling. 


Figure  3.  The  FEM  solution  for  Example  1. 


The  boundary  conditions  adopted  along  the  edges 
of  the  device  are  the  same  as  those  used  by  [16],  [12], 
and  are  displayed  in  Figure  5.  On  the  oxide- 
semiconductor  interface  the  following  boundary 
condition  is  assumed 


£ 


ox 


Vg-*-Vfb_( 


By 


(26) 
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where  Eox , and  Es  are,  respectively,  the  oxide  and 
silicon  permittivities,  VFB  is  the  flat  band  voltage,  and 
tox  is  the  oxide  thickness. 

To  solve  the  above  elliptic  problem  we  proceed 
first  by  dividing  the  boundary  into  M segments.  Finer 
segments  are  used  on  the  top  edge  of  the  device  to 
handle  its  boundary  discontinuities.  Using  the 
boundary  conditions,  and  the  technique  described  in 
[9],  matrix  [vty-]  is  computed,  then  inverted  and  stored. 
Independently,  the  base  region  is  also  divided  into  N 
finite  elements  or  cells.  We  seek  to  determine  the 
electric  potential  at  each  midpoint  of  the  cells  using 
(25)  which  is  nonlinear.  To  solve  it  we  set  up  an 
iterative  procedure,  based  on  Newton’s  linearization 
[15],  At  the  k-th  iteration  we  replace  the  right-hand 
side  of  (25)  by  its  Fourier  expansion  about  ¥*.  Then 
the  following  Helmholtz’s  equation  is  obtained 

02\j/(£+1)  02\j/(*+1) 

dx2  dy2 

F(W(k)) 

(27) 

where  F and  G are  given  by 


(*+l) 


G 


+ e 


«(*) 1 


(29) 


The  iteration  scheme  starts  by  taking  some  initial 
guess  value  for  4/(0)  so  that  (27)  can  be  solved  for  the 
first  approximation  xPf!).  Then  4/(l)  is  used  to  find  the 
second  approximation  *F(2).  The  procedure  is  repeated 

until  the  norm  II  — xy(k)  ||  js  iess  than  a 


desired  tolerance.  However,  this  iteration  scheme 
often  diverges  [13],  [15],  and  some  damping  factor 
was  found  necessary  to  improve  the  convergence  of 
the  iterative  solution.  At  the  beginning  of  the  k- th 
iteration  step,  the  following  formula  was  used  for  all 
inner  nodes 


The  accuracy  of  the  solution  obtained  by  the 
present  technique  is  demonstrated  by  comparing  it 
with  the  FDM  (5-point)  solution.  For  both  methods, 
MOM  and  FDM,  we  let  the  values  of  the  electric 
potential  to  be  updated  at  each  mesh  point  by  means  of 
an  explicit  formula,  that  is,  without  the  solution  of 
simultaneous  algebraic  equations  [15].  We  assign  *¥  = 
0,  as  an  initial  guess  for  all  the  inner  nodes.  Once  the 
convergence  is  attained  for  the  domain  nodes,  then  the 
electric  potential  at  any  other  point  in  the  device  is 
determined  by  a matrix  multiplication  as  shown  in 
(17),  and  without  the  need  of  any  interpolation.  For 
numerical  computations,  the  following  data  were  used: 
The  thickness  of  the  oxide  layer  was  tox  = 0.5  (im,  the 
flat  band  voltage  is  VFB  = -1  V and  the  doping  profile 
is  assumed  to  be  uniform  with  NA  = 1018  cm"3,  = 

1.5xl0,n  cm"3,  thermal  potential,  Vr  = 0.0258  V,  and 
R = 0.1.  The  results  of  the  computations  are  shown  in 
Figure  6,  where  the  distribution  of  the  electric 
potential  at  the  thermal  equilibrium  is  plotted  along 
different  lines  parallel  to  the  x-axis.  The  solid  lines 
represent  the  solution  obtained  using  the  FDM  and  (o) 
symbol  is  reserved  for  the  solutions  obtained  by  the 
present  method.  Close  agreements  can  be  observed 
between  the  two  methods.  However,  for  FDM  721 
nonuniform  mesh  points  are  employed  to  reduce  the 
total  number  of  nodes  {for  uniform  meshes  over  4900 
(70x70)  nodes  ought  to  be  used}.  Finer  meshes  were 
chosen  in  the  depletion  region  and  near  the  junctions 
to  reproduce  accurately  the  fast  variation  of  the 
electric  potential  [16]  and  the  surface  discontinuities 
of  the  potential;  whereas  for  MOM,  the  number  of  the 
uniform  meshes  was  529  (23x23)  or  400  for 
nonuniform  cells.  The  number  of  iterations  required 
for  the  convergence  was  5 1 with  the  present  technique 
as  compared  to  108  iterations  with  FDM  to  reach  the 
point  at  which  the  absolute  maximum  between  two 
subsequent  iterations  was  less  than  0.05  (tolerance). 

Next  we  apply  this  method  to  the  solution  of  the 
cutoff-frequencies  of  various  waveguides. 

Vo!*.}* 

(V) 


-"■’o  O 1 »).2.  VOi" CM o'i  o.ft  ' 0.7  7x» 'iy.'i j 

(Ittanc*. 


\p(fc+i)  = (}  _ z?)^^  “ o -|- 
where  R(<  1)  is  the  relaxation  factor. 


(30) 


Figure  6:  The  distribution  of  the  electric  potentials  along  the  lines 
parallel  to  the  x-axis. 
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IV.  SOLUTION  OF  THE  GENERAL 
HELMHOLTZ  EQUATION  FOR 
HOMOGENEOUSLY  FILLED 
WAVEGUIDES 

In  waveguides,  solution  of  the  Helmholtz 
equation  determines  the  electromagnetic  field 
configuration  within  the  guides.  It  is  convenient  to 
divide  the  possible  field  configurations  within  the 
waveguides  into  two  sets,  namely  TM  waves  and  TE 
waves,  each  of  which  is  governed  by  similar 
Helmholtz  equations. 

If  we  consider  a waveguide  in  which  the 
direction  of  propagation  of  the  wave  is  along  the  z- 
direction,  then  the  Helmholtz  equations  are  as  follows. 

TM.  Case  (Hz  = 0): 

V2£z  (x,  y ) + {a2  fie  - K\  )Ez  (x,  y)  = 0 (31) 

with  the  appropriate  boundary  conditions  Ez  = 0,  on 
the  conductor  walls  and  for  the 

TEZ  Case  ( Ez  s 0): 

V2Hz  (x,  y)  + (O)2 fie  - K2  )HZ  (x,  y)  = 0 (32) 

with  appropriate  boundary  conditions  dHJdn  = 0,on 
the  conductor  walls  where  dHzfdn  represents  the 
normal  derivative.  Here, 

Ez  z-component  of  the  electric  field; 

Hz  z-component  of  the  magnetic  field; 

to  angular  frequency  = 2nfi ; 

/ frequency  of  interest; 

(i  permeability  of  the  homogeneous 

medium; 

£ permittivity  of  the  homogeneous 

medium; 

Kz  propagation  constant  in  the  z- 

direction. 


Comparing  (31)  and  (32)  with  (1),  and  also 
the  boundary  conditions  of  the  TMZ  and  TEZ  cases  with 
those  of  the  general  equations,  we  can  draw  the 
following  analogies  as  outlined  in  Table  I. 

By  examining  (32)  that 

r'ti-,,  ]-[«,])(*;] 

it  can  be  inferred  that  [B]  = 0,  since  F = 0 and  y = 0 for 
TMZ  and  TEZ  cases. 

In  the  case  of  TMZ,  (22)  reduces  to  the  form 

[A]  • [Ezi]  = 0 (33) 

In  the  case  of  TE;.,  (22)  reduces  to  the  form 

[A]-  [Hzi]  = 0 (34) 

Here  again,  Ezi  and  Hzi  refer  to  the  values  of  Ez  and  Hz 
at  the  midpoints  of  the  subregions  of  the  discretized 
waveguide  cross  section. 

For  (33)  and  (34),  nontrivial  solutions  exist 
for  [E^]  and  [Hzi]  only  if  the  matrix  [A]  is  singular. 
The  condition  for  nontrivial  (i.e.,  nonzero)  solutions  to 
exist  for  [ Ezi ] and  [i/Zf]  it  is  essential  that 

det[A]  = 0 (35) 

where  det[A]  stands  for  determinant  of  [A].  We  know 
from  (23)  that 

and  we  also  know  that  for  the  cases  of  the  TMZ  and 
TEZ  waves 

X = co2  fie-K2z.  (36) 


Table  I 


Special  cases  of  the  General  Helmholtz  Equation 


General  Equation  (3) 

TMZ  Equation  (1) 

TEZ  Equation  (2) 

V2'F+?i4/=F 

V2Ez+(<n2\ui-Kl  )-Ez=  0 

V2HZ+((£>2[LE  -k\  ) Hz=  0 

Ez 

X 

w2(i  e - K2 

co2jll  e - K\ 

F 

0 

0 

a'F+pa'P/a«=7 

t»1 

(M 

II 

o 

dHz/dn= 0 

7 

0 

0 
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Hence,  given  the  frequency  at  which  the  Helmholtz 
equation  is  to  be  solved,  det  [A]  would  be  a function  of 
Kz , the  roots  of  which  will  provide  the  values  of  Kz  for 
which  det  [A]  = 0.  Once  these  Kz  values  are  known, 
the  eigenvector  of  [A]  corresponding  to  the  minimum 
eigenvalue  gives  the  nontrivial  solutions  for  [ Ezi ] and 
[Hzi]  in  case  of  TM~  and  TEZ  cases,  respectively. 

Once  [£,(]  and  [Hzi]  are  determined  at  the  grid 
points,  using  Gaussian  elimination  for  instance,  the 
values  of  Ez  and  Hz  at  any  other  point  can  be  obtained 
using  ordinary  matrix  multiplications,  as  explained  in 

tl]. 


IVa.  Calculation  of  Propagation  Constants  for 
Different  Waveguide  Modes 

It  is  evident  that  for  the  existence  of  nontrivial 
(nonzero)  solutions  for  [Ezi]  and  [Hzi],  it  is  necessary 
that  (35)  be  satisfied.  Let  us  define  a matrix  [Z]  such 
that 


[z]=([?/1][«',r'[fcjl]-[9?])  <37) 

Hence,  (35)  becomes 


Therefore,  the  propagation  constants  of  different 
modes  in  the  waveguide  are  given  by  the  following: 

For  (A:!)2 >0, 


*:=>V  e + 


For  (Ki)2< 0, 


EV 


\z  1 


K'=JJ(02ll£ 


EV}Z] 


propagating  modes 

(42) 

nonpropagating  modes 

(43) 


Results  of  propagation  constants  of  various  modes  in  a 
rectangular  waveguide  computed  by  this  method  are 
shown  next. 


IV  b.  Calculation  of  Cutoff  Frequencies  for 
Different  Waveguide  Modes 

The  cutoff  frequencies  for  the  various  propagating 
modes  in  the  waveguide  are  given  by 


det ([Z]  [A,,]  + [/])  = 0 

which  can  be  rewritten  as 


f-0 

>1 

det 

[Z]- 

[/] 

V 

J 

(38) 


(39) 


Equation  (39)  is  similar  to  the  characteristic  equation 
of  matrix  [Z],  with  its  eigenvalues  given  by 

Knowing  that  Z,  above  =(02(I  £ — K2  for  TMZ  and 
TE.  cases,  it  can  be  concluded  that 

— = 1 = RV[Z] 

A ( Klf-co 2 fie  ' ’ (40) 

i'  = l,  2,...,N. 


where  K'z  is  the  propagation  constant  of  the  ith  mode 
and  EVj^  = ith  eigenvalue  of  [Z]. 


Equation  (40)  can  be  rearranged  as 


(Ki)2 


co fie+ 


1 


EV}Z]  ' 


(41) 


/;'  = ~ /(0)2|i  £ - (Kf  ) (44) 


where  flc  = cutoff  frequency  of  the  ith  mode.  Here, 
\)  = velocity  of  light  in  the  homogeneous  medium 

= 1 / £ and  it  can  be  deduced  from  (41)  that 


CO 


AM*02=- 


-i 


EV}Z]  ’ 
t = 1, 2,...,N 


(45) 


Using  (45)  in  (44),  we  find  the  cutoff  frequencies  for 
the  first  N propagating  modes  as 


f 


-i 


EVlZ]  ’ 


i = l,2,...,N . 


(46) 

The  cutoff  wave  number  k’c  of  the  ith  mode  can  be 
calculated  from  the  cutoff  frequency  using  the  relation 

2 nfc 


k‘=- 


V 


t = l,  2,...,N 


This  method  thereby  provides  a straightforward 
approach  to  find  the  cutoff  frequencies  (and,  hence, 
cutoff  wavenumbers)  of  any  waveguide  structure 
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without  resorting  to  scanning  over  a wide  range  of 
frequencies,  as  is  done  in  the  Ritz-Galerkin  and 
surface  integral-equation  methods. 

V a.  Results:  Rectangular  Waveguide 

Consider  a rectangular  waveguide.  For  the 
waveguide  in  Figure  7,  the  region  was  divided  into 
100  subregions  and  the  boundary  was  discretized  into 
96  subcontours.  The  maximum  matrix  size  involved 
in  the  computations  was  100  x 100.  Results  have  been 
displayed  in  Table  II  for  the  cutoff  wave  numbers  of 
the  first  eight  TM/TEZ  modes.  The  computational 
time  involved  in  finding  the  cutoff  wavenumbers  of 
the  first  100  modes  on  a Sun  SPARC  10  workstation 
was  16  s. 


4 cm 

Figure  7.  A rectangular  waveguide. 


Table  II 


Vb.  Single-Ridge  Waveguide 

A single  ridge  waveguide  is  a popular  means  of 
getting  higher  bandwidth.  The  first  four  TMz  and  TEz 
mode  cutoff  wavenumbers  were  computed  for  the 
single-ridge  hollow  waveguide  shown  in  Figure  8. 
Results  have  been  displayed  in  Table  III  and  compared 
with  published  data.  For  the  waveguide,  the  region 


was  divided  into  96  subregions  and  the  boundary  was 
discretized  into  112  subcountours.  The  maximum 
matrix  size  involved  in  computations  was  96  x 96. 
The  computational  time  involved  in  finding  the  cutoff 
wavenumbers  of  the  first  96  modes  in  each  case  on  a 
Sun  SPARC  10  workstations  was  18  s. 


1.0  cm 


Figure  8.  A single-ridge  waveguide. 


VI.  CONCLUSION 

An  efficient  technique  based  on  MoM  formulation 
for  solving  a general  Helmholtz  equation  starting  from 
Laplace’s  equation  is  presented.  The  main  feature  of 
this  new  formulation  is  that  the  boundary  conditions 
are  satisfied  independent  of  the  discretizations  of  the 
regions  and  the  nodes.  This  feature  was  found 
especially  useful  when  the  boundary  conditions  have 
discontinuities.  Considerable  reduction  in  the  domain 
grids  is  realized  with  the  present  method  as  compared 
to  the  conventional  methods  such  as  finite  difference 
method  or  the  finite  element  methods.  In  addition,  one 
need  not  use  a frequency  dependent  Green’s  function 
which  can  reduce  the  computational  cost  significantly. 
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